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Abstract. In melt-blowing very thin liquid fiber jets are spun due to high-velocity air streams. 
In literature there is a clear, unsolved discrepancy between the measured and computed jet atten- 
uation. In this paper we will verify numerically that the turbulent velocity fluctuations causing 
a random aerodynamic drag on the fiber jets - that has been neglected so far - are the crucial 
effect to close this gap. For this purpose, we model the velocity fluctuations as vector Gaussian 
random fields on top of a k-e turbulence description and develop an efficient sampling procedure. 
Taking advantage of the special covariance structure the effort of the sampling is linear in the 
discretization and makes the realization possible. 
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1. Introduction 

Melt-blowing is a process for manufacturing very thin thermoplastic fibers, whose commercial 
importance steadily increases. It dates back to Wente's work in the 1950s at the Naval Research 
Laboratory in the USA [48] . For an overview on the technology we refer to (36] [28] . In a melt- 
blowing process, a molten stream of polymer is extruded from the spinneret into a forwarding high- 
velocity air stream. The aerodynamic force rapidly attenuates the polymer jet from a diameter do of 
approximately 500 micrometers at the nozzle down to final diameters d that can be as small as 0.5 
micrometers. The speeds are very high. Since air and polymer are nearly of the same temperature, 
the gas prevents polymer solidification at distances close to the die. So fibers are produced that 
are orders of magnitude smaller than the fibers of a conventional melt-spinning process where the 
stretching is caused by a mechanical force due to a take-up wheel. The elongation in melt-blowing 
is e = Ao/A = dp/d 2 ~ O(10 6 ), that means a reduction of 10 3 in diameter d and of 10 6 in 
cross-sectional area A. Melt-blown fibers make excellent filters. They have a high insulating value, 
moreover they show high cover, surface area and potentially high strength per unit weight. 

The optimization of the fabric and the manufacturing process requires the understanding of 
the fiber structure development in melt-blowing [8 . To gain insight in fiber jet attenuation and 
cooling several on-line measurements have been performed during the last years (see e.g. [511 154j 
and for measurements on jet diameter and temperature 0H7], jet velocity components [50], fre- 
quency and amplitude of jet vibrations [113 US], nonwoven webs [53] etc.). But, until now there is 
a clear, unsolved discrepancy between experiments and mathematical models / simulations. The 
numerical results presented in the literature coincide quite well with the measurements under condi- 
tions of a conventional melt-spinning process with moderate elongation e ~ O(10 2 ), but absolutely 
underestimate the jet attenuation in orders of magnitude for industrial melt-blowing processes, 
|4"Tl I5TI 155] . The reason might lie in the fact that the underlying mathematical models have been 
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Figure 1.1. Simulation of turbulent spinning in a melt-blowing process. Ver- 
tically downwards directed high-velocity air stream with immersed spun fiber jets 
(white curves) whose random motion is caused by the turbulent fluctuations. The 
color scale visualizes the vertical component of the mean flow velocity. For details 
see Section |U 



originally developed for melt-spinning processes, dealing with mass, momentum and energy balances 
for a steady (longitudinal) spinning threadline, cf. first publications [T5J [33J in the 1960s or for an 
overview 56J. Up to now the studies have been extended to viscous and viscoelastic fluids with 
inclusion of heat transfer, inertial and air drag effects and with regard of jet dynamics, vibrations 
and bending instabilities. It is an area of active research as recent articles show, see for example 
[T3], E21 EOl EU EH1 E3l [44J, [52l E3] and references within. However, in the used steady considerations, 
the jet cross-sectional area A and the speed v are related according to voAo = vA for an incompress- 
ible fluid where the index indicates the quantities at the nozzle, such that the computed elongation 
is obviously restricted by the velocity u of the acting air stream, i.e. e = v/vo < HuUoo/vo. This 
also holds true for advanced melt-blowing simulations with a turbulence model for the high-velocity 
air stream when only the mean flow field informations are taken into account |55j . The computed 
elongation is hence e ~ O(10 4 ) - in contrast to the measured e ~ O(10 6 ). Latest experiments 
[5T1 S3J indicate the relevance of the turbulent fluctuations for the jet thinning. 

This paper aims at the numerical verification of the crucial effect of the turbulent velocity fluctu- 
ations for the fiber jet attenuation (Figure [L~T| ) . The turbulent fluctuations of the high- velocity air 
stream turn out to be qualitatively very important, their impact might close the gap between mea- 
surements and previous simulations that has been observed and studied for a long time. Since direct 
numerical simulations (DNS) are computationally not possible, we use a stochastic k-e turbulence 
description for the high-velocity air stream, on top we model the turbulent fluctuations as Gauss- 
ian random fields in space and time according to [291 131) . Their covariance/correlation functions 
satisfy Kolmogorov's 5/3-law for the energy spectrum and the requirements of the k-e turbulence 
model. Moreover, they have a special structure, for which we develop an efficient sampling strategy 
that makes the application possible. Under the conditions of an industrial melt-blowing process we 
demonstrate the relevance of the turbulent fluctuations and the resulting random aerodynamic drag 
force in a one-way-coupling between the air stream and the fiber jet. Thereby, we use an isothermal 
model for the fiber jet dynamics that consists of a system of ordinary differential equations (ODE) 
for jet position, velocity and elongation. Already for this very simplified ODE-model with random- 
ness we obtain a jet thinning behavior that is qualitatively appropriate in magnitude. The results 
are very promising and raise hope that the application of the random aerodynamic drag to more 
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sophisticated Cosserat models for the jet dynamics [H 21 [3pJ SO] that also include inner stresses and 
temperature dependencies (systems of partial differential equations (PDE)) will finally answer the 
open questions of the fiber structure development in melt-blowing in future. Recent work deals with 
the robust numerical treatment of the PDE- models [3J |S] . Of particular importance is thereby the 
establishment of a fast and accurate adaptive mesh refinement that is necessary in order to cope 
with the expected huge elongations. 

The paper is structured as follows. After a short survey in turbulence modeling, the reconstruc- 
tion of the turbulent velocity fluctuations as Gaussian random fields taken from [29 is presented in 
Section[2] In Section|3]we develop a suitable, very efficient sampling strategy that makes use of the 
special covariance structure. We apply the resulting random aerodynamic drag force to a simplified 
ODE-model for melt-blowing in Section [4] Comparing the results with regard and neglect of the 
turbulent fluctuations clearly shows the importance of the fluctuations for a proper description of 
the manufacturing process. 

2. Turbulence Modeling 

All the physics of a turbulent flow is contained in the non-stationary Navier-Stokes equations 
(NSE). Turbulence can be considered as continuum phenomenon, since even the smallest scales 
occurring in a turbulent flow are ordinarily far larger than any molecular scales. Thus, solving NSE 
by means of DNS gives the exact velocity field needed for the determination of the aerodynamic 
force causing the jet stretching. However, DNS presupposes the resolution of all vortices ranging 
from the large energy-bearing ones of length It to the smallest, viscously determined Kolmogorov 
vortices of size ij with \t/t) — Re 3//4 [15] ■ Hence, the number of grid points that are required 
for the refinement of a three-dimensional domain is proportional to Re 9 ^ 4 . Despite of recent high 
speed computations, DNS is thus still restricted to simple, small Reynolds number flow. In large 
eddy simulations (LES) the computational effort is reduced. By applying a low-pass filter on NSE, 
only the vortices of large scales are resolved directly, whereas the small vortices are taken into 
account by a stochastic approximation of their effect on the larger ones. This procedure works well 
for fluid-structure interactions with structures of moderate size, e.g. flow around the wing of an 
aircraft. However, in view of our thin fiber jets we need correlation information and characteristics 
of such small scales so that LES is not applicable. Therefore, we use stochastic turbulence models. 
These models suffer from their deficient description of boundary layers at walls and/or obstacles. 
But since we restrict to a one-way-coupling by neglecting the effect of the few slender jets on the 
air flow, the flow information that we need for the air drag model comes from the simulation of 
a turbulent flow in a machine geometry without any immersed fiber jets. So the vortex structure 
develops in the inner flow domain, unperturbed by any obstacles. And no inaccurate boundary 
layer approximations falsify the numerical results. The stochastic turbulence models - also known 
as statistical turbulence models [HI 03] ~ prescribe the instantaneous space- and time-dependent 
velocity field u : K 3 x K^" — > K 3 as sum of a mean (deterministic) u and a fluctuating (stochastic) 
part u' 

u = u + u'. 

So, the velocity field is considered as a R 3 -valued random field u = (u(x, t))^ x t j eR 3 xR +, i.e., the 

velocity at the point (x, t) in space and time is assumed to be a realization (sample) of the Un- 
valued random variable (random vector) u(x, t). On the average over all realizations u equals u, 
i.e., E(u) = u or, equivalently E(u') = 0, where E denotes the expectation. The mean velocity 
u is determined by the Reynolds-averaged Navier-Stokes equations (RANS) where the effect of 
the fluctuations is contained as source in terms of the Reynolds stress tensor T = — E(u' eg u') 
in the momentum balance. The fluctuations u' = (u'(x, t))^ x t ) eR 3 xR + are not modeled directly 
themselves as centered random field, instead the stochastic turbulence models restrict on modeling 
the symmetric Reynolds stress tensor. Therefore additional statistical quantities are introduced. 
In the k-e model, one of the prime turbulence models going back to Launder and Sharma [23], 
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these quantities are the turbulent kinetic energy k : M 3 x Kg — » R + and the dissipation rate 
e : M 3 x Rq — > M+. They characterize the fluctuations 

k= -E(u'-u'), e = i/E(Vu':Vu') (2.1) 

with flow viscosity z/ and yield the Reynolds stress tensor T = C M k 2 /e (Vu + Vu T ) — 2k/3J with 
unit tensor I and constant C M in accordance with the Boussinesq assumption that proceeds from 
the analogy between viscous and turbulent stresses. In general, the purely deterministic PDE- 
models for the turbulence description consist of RANS and transport equations for the additional 
statistical quantities. They are closed by means of numerous assumptions and fitted parameters 
(closure relations). 

Notation 1 (Tensor Calculus). A tensor can be associated with a multi- dimensional array with 
regard to the chosen basis, e.g. a tensor of 2nd, 1st or Oth order is represented by a matrix, vector 
or scalar, respectively. Throughout the paper we use large and small bold-faced letters for matrix- and 
vector-valued quantities, respectively. Scalar-valued quantities are typeset in normal-faced letters. 
Elementary operations include tensor and dot products, such as a ® b = C with Cy = a^bj as well 
as a • b = &ihi and A : B = ■ Ay By (scalar products of vectors and matrices). 

On top of a k-e formulation we developed a turbulence reconstruction strategy for u' in 29J . This 
strategy is based on a Global-from-Local- Assumption according to that the local velocity fluctua- 
tions (fine-scale structure) are modeled as homogenous, isotropic Gaussian random fields that are 
superposed to form the large-scale structure of the global turbulence. This assumption is motivated 
by Kolmogorov's local isotropy hypothesis |15j : certain theoretical considerations concerning the 
energy transfer through the eddy-size spectrum from the larger to the smaller eddies lead to the 
conclusion that the fine structure of anisotropic turbulent flows is almost isotropic. 

Assumption 2 (Global-from-Local- Assumption [29). Let the k-e turbulence model be given. Let 
(Cl,A, P) be a probability space. Let u'; ocp = (u'; oc . p (x, t))^ x t ) eR 3 xR + be a family of centered Gauss- 
ian random fields on (fi, A, P) that correspond to spatially and temporally homogeneous, isotropic, 
and incompressible flow fluctuations with respect to the local turbulence information (k, e, v, u) at 
the point p = (y, t) € R 3 x Rq~. Their tensor-valued covariance/ correlation functions are denoted 
by K/ OCi p : (R 3 x Rq ) 2 — > R 3x3 . 

Then we assume that the actual global fluctuation field u' can be constructed as 

u'(x,t) = (u' /oc , p (x,t)) M(x , t) , (2.2) 

withM{x,t) = {p= (y,r) el 3 x R+ | ||x - y - u(x, t)(t - r)|| 2 < 1 T A |t - r| < t T } ; |M(x,t)| = 
/m(x t) ^P' an< ^ turbulent large-scale length It and time tx- The brackets (•) stand for a Gaussian 
average with respect to the parameter^, i.e., u' is a Gaussian random field that is uniquely prescribed 
by the following mean and covariance function 

E(u'(x, t)) = / / E(u' ioc , p (x, t)) dp = 0, 

\ M \ X ^)\ JA/(x,t) 

E(u ; (x,t) O u^xi^i)) = 1 == / K Joc , p (x,t,xi,ti)dp 

A/|M(x,t)||M(xi, tx) | JM(x,t)nM(xi,ti) 

= K(x,t,x 1 ,t 1 ). 

Remark 3. A random field is Gaussian if all its finite- dimensional joint distributions are Gaussian 
(normal distributed). The distribution of a Gaussian random field is completely specified by its mean 
and covariance function. By Assumption^ the modeling of the turbulent fluctuations u' is reduced 
to the modeling of u[ oc focusing on an appropriate description of the covariance function K; oc p . 
The covariance function contains the information about the spatial and temporal correlations in the 
turbulent flow. 
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The local centered Gaussian velocity fluctuation fields u'; OC)P and hence their covariance/correla- 
tions K; oc p depend parametrically on the flow situation (k, e, v, u) at the point p that is provided 
by the k-e model. Thus, we make u'; OCiP dimensionless using the typical turbulent length k 3 / 2 /e 
and time k/e corresponding to p, 

/ / \ ,1/2/ / 6 e e \ k 3 / 2 k k 2 , 

u loc (x,t)=k ' n ioc ^j-^x, -t; ^i/J , x=— — x, t = -t, v = — f. (2.3) 

The dimensionless viscosity £ enters the model of the velocity fluctuations via the consistency with 
the k-e description, see Model [5] To facilitate the readability we suppress the parameter-dependence 



(index p) here in (2.3) and also in the following. 



Notation 4 (Dimensional vs Dimensionless Quantity) . We typeset dimensional quantities in Roman 
style (e.g. t, x, u' ; K.) and the corresponding dimensionless quantities in Italic style (e.g. t, x, u' , 
K ) throughout the paper. 

In the forthcoming explanations we focus on the dimensionless quantities. The development of 
the local tensor-valued covariance/correlations K i oc can be reduced to the modeling of two scalar- 
valued functions, the energy spectrum E and the decay of the temporal correlations ip. (For the 
detailed derivation of the original model with frozen turbulence pattern we refer to |29j , extensions 
on the temporal correlations are studied and incorporated in [31 j - For more informations about 
turbulence and its evolution see [TH1 [T71 [27J [3T] and references within.) In a homogeneous turbulent 
flow, the correlations are invariant with regard to spatial and temporal translations and hence 
depend only on the differences of the arguments. The evolution of the correlations are modeled 
by an advection-driven vortex structure that is naturally decaying over time (alleviated frozen 
turbulence). Taylor's hypothesis of frozen turbulence pattern [45], i.e. fluctuations arise due to so- 
called turbulence pattern that are transported by the mean flow without changing their structure, 
is based on the observation that the rate of decay of the mean properties is rather slow with 
respect to the time scale of the fluctuating fine-scale structures. The superposition with a natural 
temporal decay is essential for describing suspensions of particles or filaments in turbulent flows, 
since otherwise small light objects tending to move with the mean flow field would experience 
permanently the same non-varying fluctuations. So, K i oc is prescribed by the initial correlation 
tensor 7 : R 3 — > R 3x3 and the temporal decay function <p : R+ -> R, 

K loc (x + x 1 ,t + t 1 ,x 1 ,t 1 ) = j(x - ut) ip(t) (2.4) 

with mean flow velocity u. It is also made dimensionless with respect to the typical turbulent length 



and time in consistency to (2.3 1. In case of incompressible isotropic turbulence, 7 is represented by 
a single scalar- valued smooth function. In terms of the spectral density being its Fourier transform 
J 7 ^, this function is known as energy spectrum E : R$ — > R^ that has been well-studied theoretically 
and experimentally in the last century, 

1 ^(IMI) (r 1 . « 



F~({k) = [ exp(— in ■ x) -y(x) dx 



4tt 



(2.5) 



with unit tensor J and Euclidian norm ||.||. Kolmogorov's universal equilibrium theory was thereby 
trend setting. Based on dimensional analysis he derived not only the characteristic ranges but 
also the typical behavior of the spectrum which agrees with later coming physical concepts and 
experiments, cf. Kolmogorov's 5/3-law and his hypothesis of local isotropy [TS]. Gathering the 
existing knowledge about E, an appropriate model has to satisfy Kolmogorov's 5/3-law as well as 
the requirements of the k-e turbulence model, i.e. 

/ E(k) d«=l, / K 2 E(n)dn=— , / (ln(l + k)) q k 2 .B(k) &k < 00 for some a > 3. 
Jo Jo 2C Jo 

The first two relations correspond to the definitions of the kinetic turbulent energy k and dissipation 



rate e in (2.1 1. For e to make sense, the third relation ensures that the fluctuation field is almost 



surely sample diffcrcntiable in space. (The n-times sample differentiability of a Gaussian field 
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results from certain integral properties of its spectral (density) function, for details see [3T].) The 
conditions on E allow for a family of functions that can be adapted to experiments. The parametric 
C-dependence comes from non-dimensionalizing (2.1) and is handed over to the correlation tensor 
and u'loc, cf. ( [273] ). 

Model 5 (Energy Spectrum [29l El] ) - The energy spectrum is modeled as E 6 C 2 (lRg~) 

H 5/3 E-= 4 Orf & « < «1 
C) = C K { k" 5/3 ki < k < k 2 (2.6) 

where the (^-dependent transition wave numbers ki and K2 are implicitly given by 

/ E(«;OdK=l, / K 2 £( K ;C)d« = — . 

The regularity parameters are = 230/9, 05 = —391/9, a§ = 170/9, 67 = 209/9, 6g = —352/9, 
bg = 152/9, and the Kolmogorov constant is Ck = 1/2. 



The integral conditions for K\ and k 2 in £ (2.6) can be reformulated as nonlinear system 

« -2/3 f -2/3 ^-1 . 4/3 . f 4/3 /nr , /„ „S 

ai - 2 + ^JTT' a2 ~4~^JT3' fel - 2 "^J3T' 62 - 4 + ^J^3- 

J= 4 J j=4 J j=7 J j=7 J 

The condition < K\ < n 2 < 00 is equivalent to < ( < (crit = (2C|-(&2 — 02X^1 — Si) 2 ) -1 s» 3.86. 
The bounds on £ (where we have n\ — n 2 — (Ck{cli — fri)) 3//2 for ( = Q r a and ki = (Ckcii) 3 ^ 2 , 
K2 = 00 for £ = 0) are no practically relevant restrictions, since the general turbulence theory 
assumes the ratio of fine-scale and large-scale length to satisfy ( = ev/k 2 <C 1. 

The temporal correlation function ip satisfies ip(0) = 1 which implies that the integral of its 
Fourier transform T v is normalized. We use an exponential decay with respect to the turbulent 
large-scale time with tj- = 0.212, see e.g. [25J|57] and references within. 

Model 6 (Temporal Correlations). The natural decay of the temporal correlations is modeled as 

t^ f — t^uP 1 



Remark 7. The correlated global random field u' (2.2) can be asymptotically reduced to Gaussian 
white noise with flow- dependent amplitude. For certain applications this simplification is quali- 
tatively and quantitatively justified and implies an enormous reduction of computational time and 
memory ( see |29j for a theoretical localization strategy and the strict asymptotic derivation as well as 
|31j for the application and experimental validation in a production process of nonwoven materials). 



3. Sampling of Gaussian Random Fields 

In this section we deal with the reconstruction and simulation of the local centered Gaussian ran- 
dom velocity fields u' ' \ oc . In view of the prescribed data there exist various reconstruction/sampling 
procedures in literature, e.g. Karhunen-Loeve expansion, Cholesky decomposition, circulant embed- 
ding for a given covariance function [9| or spectral, Fourier, Fourier-wavelet methods for a given 
spectral function [TT] EE H2] ; f° r an overview see [5J [2U] and references within. We propose a tech- 
nique that takes advantage of the special structure of the given data and turns out to be exact in 
the covariance and very efficient as we will comment on. Note that in the following we assume the 
existence of all occurring random fields and stochastic processes as we construct them later on. 
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3.1. Construction of u'i oc . Considering the covariance function (2.4) of the local velocity field, 
we separate the spatial and temporal arguments by introducing a new Gaussian random field r/ = 
(»70M))(x,i)en 4 

T](x,t) = ll' loc (x + Ut,t) 

from which u'i oc can be easily regained. Its covariance satisfies 

E(rj(x + x 1 ,t + t 1 )® riixt , h)) = 7(35) <p(t). (3.9) 

Let the vector random field £ = (£(x)) xS R3 and the scalar stochastic process ip = (ip(t))teM be 
centered and stochastically independent with covariance functions 

E(£(x + x J )®£(x i ))= 7 (x), E(V>(t + *i)V(*i)) =¥>(*)• 

Defining a random field i) by 



i) and r\ possess the same covariance function ( |3.9| ). As we are interested in a Gaussian field, we 
consider 

1 N 

f lN (x,t) = ^= p ^r l «\x,t), NeN, (3.10) 



in which f)^ l \ I — 1, ...,N are independent copies of r). The central limit theorem ensures then the 
convergence in distribution 

for every (x, t) £ K 4 as N tends to infinity. Due to the multi-dimensional central limit theorem, for 
any choice ofneN and (xj , ii), . . . , (x n ,t n ) e R 4 , the joint distribution of fi N (x\,t\), . . . , f] N {x n ,t n ) 
converges in distribution to a normal distribution on M. 3n . We conclude that f/ N is a centered ran- 



dom field with covariance (3.9), which is approximately Gaussian if N is large. So in order to 
construct tj respectively f} N we focus on the sampling of £ and ip. Thereby, we keep in mind the 
almost sure differentiability of the realizations of our constructed fields and processes. 

3.1.1. Spatial Field In this subsection we exploit the special structure of the spectral density of 



the spatial field given by = in (2.5). Let w = (u>(i)) 4eH be a centered, homogeneous and 
Revalued stochastic process with spectral function S w (k) = E(\k\)/2I, i.e., its components are 
uncorrelated processes with the same spectral function s w (k) — E(\k\)/2. Moreover, let z be a 
uniformly distributed random vector on the unit sphere S 2 = {x € E 3 : ||x|| = 1}. Then, under the 
assumption that w and z are independent, the random field that is defined by 

£(x) = (Z — z ® z) ■ w(x ■ z) 



has the spectral density S% given by ( |2.5[ ) and hence the desired covariance 7, [261 [12]. Since the 
components of w are independent, it is sufficient to focus on the sampling of one component w in 
order to construct the whole field Following [22 , this can be done in the subsequent manner: 
As E(k) > for all k > and J ffi s w (k) d« = J Q E(k)&k, = 1, the function s w is a continuous 
probability density on K. Choosing a random variable R with this probability density and two 
standard normally distributed random variables X and Y that are all stochastically independent, 
the complex-valued process (w(t))teR 

w(t) = Zexp(iRt), Z = X + iY, 
has the spectral function 2s w as a simple calculation with the conjugate-complex w shows 

E(w(t + ti) w(h)) = E (exp(iRt)) E(ZZ) =2 exp(iKi) s w (k) dn. 
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By taking its real or imaginary part we obtain a real-valued process with the desired spectral 
function s w . The so constructed process w (w = Re(w) or w — Im(to)) has obviously almost surely 
diffcrcntiablc realizations and hence the same holds for £. 

3.1.2. Time Process ip. The sampling of the time process ip can be performed analogously to w. 
We introduce the new process ip(t) = ip(t^t) having the covariance E(ip(t + ti)ip(t)) = (p(trt), cf. 



(|2.8|). Consequently, its spectral function is 

s ^ )= ^(£) = vt exp H 2 

As s~r is the probability density of the standard normal distribution, we take three independent, 
standard normally distributed random variables R, X, Y and set ip(t) — Z exp(iRt) with Z = X+iY. 
Then, the process ip — Ke(ip(-/ti)) or ip — Im(ip(- ft?)) has the desired covariance function tp and 
almost surely differentiable realizations. 

The reconstruction strategy for £ uses the isotropic form of the spectral density and traces the 
construction of a random field back to the sampling of a scalar-valued stochastic process which 
involves an enormous reduction of complexity and effort. For the sampling of the scalar- valued 
process w with respect to the spectral density, various (approximate) Fourier methods could be 
applied. However, our approach of introducing the complex-valued surrogate process via three 
random variables is not only exact but also efficient. The sampling and evaluation can be performed 
flexibly regarding the needs. In contrast to the a priori fixed discretization in the Fourier methods, 
this adaptivity yields advantages concerning memory and computation costs for the forthcoming 
simulations of the jets dynamics. The same holds true for the time process ip. Here, one could 
certainly think of direct methods on top of the covariance, but their performance suffers from an a 
priori discretization and high effort (for example the effort for a Cholesky decomposition is 0(n 3 ), 
n number of grid points). Our effort is linear in the discretization. For every u'i oc we have to 
generate 9N standard normally and 3iV s^-distributed variables as well as N uniformly distributed 
vectors on S 2 . Thereby, the realization of the s^-distributed variables is the most expensive part. 
These variables depend on the considered flow situation as s w (k;Q — with £ = ev/k 2 

at point p. 

Remark 8. As for the stochastic simulation (cf. Algorithm^, a random vector z that is uniformly 
distributed on the unit sphere S 2 can be sampled by help of three independent standard normally 
distributed scalars X\, X2, X3 according to z — {X\, X%, X%)/R with R = \J X 2 + X\ + X 2 (scaling 
method). For the generation of a s w -distributed variable we use the classical acceptance-rejection 
method by von Neumann |35j . taking the gamma distribution as reference density. For details on 
the techniques see e.g. [191 134] . 

Algorithm 9 (Sampling Procedure). 

Output: approximate sample from u'i oc at (x,t) 

Input: flow data at point p: k, e, v, u and dimensionless turbulent large-scale time t^, 
evaluation point (x,t) 

A.l Determine the dimensionless local flow parameters: £ = ev/k 2 and u = u/vk 

A. 2 Set random field parameter N and generate random numbers, 1 = 1,..., N : 

► z(0 

N samples according to the uniform distribution on the sphere S 2 by scaling method 

► Jig, j = 1,2,3 

3^ samples according to the density s w by von Neumann's method 
(obtain s w (n:X) — -^(l K hC)/2 by solving the nonlinear system ( |2.7[ ) ) 

► Xg, Y^, j = 1,2,3 as well as R® , X%\ Y® 

9N samples according to the standard normal distribution 

B. l Compute approximate samples, I = 1, . . . , N: 
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► spatial field: 

wf{x ■ z {l) ) = Re((xg + iY®) exp(ii?g x ■ z^)), j = 1, 2, 3 



► iime process: 
B.2 Approximate 



^ (t) = Re((X« + zT«) exp(»JiWt/*r)) 



u' loc (x, i;C)«^=E £ W ~ u*) (*) 



Algorithm [9] consists of two parts, A - the initialization with the generation of random numbers 
and B - the computation and evaluation. Hence, to evaluate the same sample at a different collection 
of points (xi,ti), only part B need to be executed while the initialization with the random numbers 
of part A should be stored. 

3.2. Simulation of u'i oc and Tests on Multivariate Normality. Simulating the local velocity 
fluctuations, the finite-dimensional distributions of fj N are close to a multivariate normal distri- 
bution for large N according to the central limit theorem. For the assessment of the multivariate 
normality on a fixed set of points {(xi,ti), . . . , (xd,td)} C M 3 x K^" we apply the statistical test 



by Royston [IT], [32] which we evaluate by help of the respective MATLAB routine [46] . Table 3.1 
shows the rejection frequencies at the significance level 0.05 for different values of the variate size d 
and the random field parameter N. We use here 1000 Monte Carlo replications and a sample size 
of 50. The rejection frequency among the 1000 replications turns out to be a robust quantity that 



N\d 


1 


2 


3 


4 


5 


6 


10 


0.357 


0.416 


0.439 


0.461 


0.497 


0.518 


30 


0.108 


0.15 


0.184 


0.187 


0.19 


0.182 


50 


0.098 


0.102 


0.12 


0.143 


0.183 


0.126 


70 


0.069 


0.079 


0.072 


0.087 


0.117 


0.124 


100 


0.082 


0.076 


0.096 


0.091 


0.085 


0.108 


150 


0.094 


0.099 


0.093 


0.101 


0.109 


0.125 



Table 3.1. Rejection frequencies of Royston's test. 




y x y 



Figure 3.2. Realization of a component of u'i oc with ( = 0, plotted over two- 
dimensional space at certain times. 
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is widely independent of the variate size d, moreover it stays approximately the same for N > 50. 



Hence, we use N — 50 in (3.101 for the forthcoming simulations. The observed rejection frequency of 



10% is acceptable for us since the Gaussian distribution is a secondary property. Our main interest 



is the accurate construction of the covariance structure. Figure 3.2 illustrates the numerical result 
of a realization of u'[ oc for a typical flow situation with £ = 0. 



3.3. Globalization Strategy Realization of the Global Velocity Fluctuations. According 
to the Global- from-Local- Assumption (Assumption [2]) , the global velocity fluctuations u'(x, t) can 
be obtained from the Gaussian average over the local information of the helds u'; ocp that belong 
to the neighborhood p £ M(x, t) of relevant space and time correlations (2.2). This construction 
satisfies the requirements of the k-e turbulence model (2.1) in an integral sense, It becomes 



already very memory-demanding and time-consuming for a discretization with |M(x, t)| > 5 which 
is obviously rather coarse in R 4 . In view of a turbulent spinning process it is not applicable. 
Thus, we propose the following globalization strategy 



u'(x, t) = k x / 2 (p) «V P (^j (p) x, ~(p) t; ~(p)) 



p=(x,t) 



(3.11) 



using the space- and time-dependent flow fields for kinetic turbulent energy, dissipation rate and 
viscosity known from the k-e-simulation. The global velocity fluctuation field of ( 3.11[ ) fulfills 
the condition (2.1) on the kinetic turbulent energy exactly, i.e. E(u'(x, t) • u'(x, t)) = 2k(x,t) for 
all (x, t) £ K 3 x Kq". Furthermore, the condition on the dissipation rate containing the spatial 
derivatives is valid up to an error of order O(Co) with constant £ = eg^o/k 2 , <C 1 representing the 
typical ratio of turbulent fine-scale and large-scale length. This estimate is based on the observation 
that changes in the behavior of k and e mainly appear on the large scale and not on the fine scale. 
This motivates an asymptotic consideration with the multi-scale ansatz k(x, t) = ko + ki(£ox, t) 
(and analogously for e, v), yielding the result. 

Since the whole field f (x, t) = ev/k 2 (x, t) is negligibly small in turbulent flows (see e.g. Figure 
for a turbulent air stream in a melt-blowing process), (3.11 1 might be further simplified to 



4.5 



u'(x,t) =k 1 / 2 (p) u' 



Ik 3 / 2 



(p)x,-(p)t;0) 



p=(x,t) 



(3.12) 



Considering an asymptotic expansion in Co, (3.12) and (3.11) obviously agree in leading order 
By this slight modification 

;C 



the sampling of the global velocity fluctuations u' in (|3.12|) can be 
performed with respect to a parameter-free energy spectrum E( ■ 



0) (2.6), which avoids the 
expensive solving of different nonlinear systems (2.1) and the multiple application of von Neumann's 
method for all the required u'i oc p (see Algorithm |9j Step A. 2). So, in total we only need a single 
set of random parameters (3N s^-distributed variables with s w (k; 0) = B(|k|; 0)/2, k £ M. as well as 
9N ~ Af(0, 1) and TV uniformly distributed vectors on S 2 ). This involves an enormous decrease of 
computational costs and makes the globalization strategy applicable to practically relevant problems 
as we will show. Note that its realization is straightforward based on Algorithm [9] 



4. Application to a Turbulent Spinning Process 

4.1. Simplified ODE-Model. The characteristics of a turbulent spinning process like melt-blowing 
are the huge jet elongations (jet thinning) that are obtained by the stretching due to turbulent air 
flows. Up to now, the numerical simulations available in the literature, e.g. [471 1551 151 ) . cannot 
predict the large elongations measured in the experiments e ~ O(10 6 ). We suppose that the reason 
lies in the steady considerations of the turbulent flow field and the neglect of the fluctuations. In 
|43j perturbations (bending instability) on a jet were imposed by turbulent pulsations, leading to 
stretching and thinning. In [2^1 OH] we developed a model framework for the dynamics of a long 
slender object (fiber) in a turbulent flow in terms of a random aerodynamic drag force in a one- 
way-coupling. The dimensionless drag force / depends on the relative velocity between air flow 
and object as well as on the object's tangent. To study the impact of our random force model and 
to get first qualitative estimates for the jet attenuation, we deal here with a simplified model for 
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an isothermal fiber jet dynamics. Instead of the complex PDE-based Cosserat models that contain 
inner stresses and temperature dependencies [TJ[21[I]j it just consists of a system of first order ODEs 
in time for jet position r, velocity v and elongation e. For this purpose, we approximate the jet 
tangent (space derivative) in the drag force / by the direction of the jet velocity v/||v|| and moti- 
vate an evolution equation for the elongation from the steady situation where e = ||v||/vo with exit 
velocity magnitude at the spinning nozzle vo (cf. Section [lj. The resulting model (4.13) describes 
the path and behavior of a single jet point whose motion is exclusively driven by a turbulent air 
flow. Certainly, gravitational forces could be included, but they play a negligibly small role for the 
observed jet thinning. The constructed Gaussian field u = u + u' for the turbulent flow velocity 
carries the randomness into the dynamic system via the drag force, 
d 



dt 
d 
dt 
d 

dt' 



v = e 3/2 a(r,t)/ 



I 



,t) 



.3/2 



V 



a(r,t) 



/ 



b(r,t) 
1 u(r,t 



b(r,t) 



r(0): 
v(0) 

e(0): 



VoT 



(4.13) 



where the scalar-valued functions a and b 

4 

a : 



pV 



b = 



V 



7T ppdg ' 

contain - apart from constant fiber jet quantities (pp density, d diameter at the nozzle) 
space- and time-dependent air flow information (p density, v viscosity). 



also 



Remark 10 (Numerical Treatment). For the forthcoming numerical investigations of the random 
ordinary differential system (4.13) the k-e- simulations of the underlying turbulent flow field are 
performed with the software ANSYS Fluent. The fiber jet dynamics is computed in MATLAB 
using the standard ODE-solver odeJ^.5.m. The routine is an implementation of an explicit Runge- 
Kutta method of fourth ( or respectively, fifth ) order with adaptive time step control. The random 
normally and uniformly distributed numbers that are required for the sampling of the turbulent 
velocity fluctuations u' are generated with the MATLAB-functions randnO and rand(). 

4.2. Numerical Results. We investigate the dynamics and behavior of a fiber jet in a flow sit- 
uation that is usual for melt-blowing, [351 ISS]- Temperature effects are neglected for simplicity. 
The air stream is directed vertically downwards and enters the domain of interest via a thin slot 
die, cf. Figure 4.3 Since the mean quantities of the turbulent air stream are time-independent 
and homogenous in direction of the slot (x-dircction), we perform stationary k-e-simulations for 
a representative two-dimensional cut showing the y-z-plane. In the set-up the mean flow is sym- 
metric with respect to the z-axis (y = 0). Figure |4~4| shows the respective flow fields for the mean 
velocity components, kinetic turbulent energy and dissipation rate; in addition p«l [kg/m 3 ] and 
v = 1.5 • 10~ 5 [m 2 /s] is constant. A distinct free air jet develops that is supersonic at the inlet 
slot (here: ||u|| 400 [m/s], k w 10 3 [m 2 /s 2 ], e « 10 s [m 2 /s 3 ]) and becomes subsonic within some 
centimeters away. So, the occurring typical turbulent length and time scales lie in a wide range, i.e., 
1 T = k 3 / 2 /e <= (1(T 4 , 10~ 2 ) [m] and t T = k/e £ (1CT 5 , 10~ 3 ) [s]. But, C ~ 10" 4 in the whole free air 



stream, as visualized in Figure |4.5| At the boundaries of the flow domain we observe side effects 
coming from the geometry (e.g. in the lower corners). However, these play no role for the dynamics 
of the fiber jet. Consequently, the simplified globalization strategy (asymptotic limit ( = (3.12)) 
for the sampling of u' is acceptable and applied. 
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Figure 4.3. Sketch of flow domain with immersed fiber jets. A two-dimensional 
cut (y-z-plane, marked by dashed line) is representative due to the given homo- 
geneity in x-direction. For details on possible geometries see e.g. [251 155] , 




Figure 4.4. k-e-simulation of the representative 2d flow domain. Top: compo- 
nents of mean velocity u in y- and z-direction. Bottom: turbulent kinetic energy k 
and dissipation rate e in logarithmic plots. 
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Figure 4.5. Turbulent scales corresponding to Figure |4.4| in logarithmic plots. 
Top: turbulent large-scale length It = k 3 / 2 /e and time tx = k/e. Bottom: ratio of 
fine and large scales C = e^/k 2 . 



The immersed fiber jet to be spun in (negative) z-direction is initialized at the slot die (spinning 
nozzle) with vo = 10~ 2 [m/s], do = 4 • 10 -4 [m] and — 7 ■ 10 2 [kg/m 3 ] and simulated for the 
time interval [0, T), T = 1CP 3 [s]. The fiber jet moves exclusively in the distinct region of the free 
air stream, Figure |4.6| Thereby, its motion is determined strongly by the mean flow pulling the 
fiber jet straight downwards, on the one hand. On the other hand the flow fluctuations cause a 
slight bouncing. The fiber velocity follows and finally adjusts to the flow velocity, as we can see 



in Figure 4.7 In the temporal evolution the fiber point starts from the nozzle where the impact 
of the turbulence is at the strongest. The velocity fluctuations act here on tiny length and time 
scales causing a quick acceleration and a very strong stretching of the jet. When the fiber point 
is some centimeters away from the nozzle after 0.2-0.3 milliseconds, the turbulence attenuates and 



the turbulent scales become larger (Figures 4.5 and 4.9). In particular, tx and the jet's reaction 
time coincide which can be concluded from the velocity curves that match. Also the elongation 
stagnates. 
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Figure 4.6. Fiber dynamics driven by turbulent aerodynamic drag force. Left: 
random trajectory r. Right: projection of several trajectories into y-z-plane (white 
curves). They are located in the distinct free air stream where C ~ 10 -4 <C 1. 
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Figure 4.8. Elongation. Top left: (sample) fiber elongation e over time with 
regard and neglect of turbulent velocity fluctuations u'. Other plots a)-c): proba- 
bility density of e at depicted heights r3 = —0.033,-0.066,-0.1 [m] (z-direction) 
estimated with Monte-Carlo simulation. 



Figure 4.8 shows the probability density function of the elongation at three depicted heights 
i\3 = —0.033, —0.066, —0.1 [m]. It is estimated by help of Monte-Carlo simulations with 5000 
samples. We clearly observe the essential effect of the turbulent velocity fluctuations and our random 
aerodynamic drag force model on the jet thinning (especially in the first centimeters / tenths of 
milliseconds). The computed elongation rises up to a mean of 2- 10 5 at T (here: r3 w —0.14 [m]). In 
comparison, the numerical result neglecting the fluctuations is approximately 10 4 which perfectly 
corresponds to the theoretical considerations on stationary turbulence stating that e = ||v|j/vo with 



v s» u holds, see Figure 4.8 The fact that Zeng et al. [55] have obtained the same magnitude of 
elongation for an visco-elastic spring-beam jet model in a mean turbulent flow field clearly stresses 
that the turbulent fluctuations are the major dominant effect for the large jet attenuation; material 
models and inner stresses in contrast seem to be of minor relevance. Moreover, it is worth to 
mention, that our ODE-model - as simple as it is - already predicts qualitatively appropriately 
all jet thinning stages observed in the experiments. However, proper quantitative estimates can be 
only expected according to the measurements [7] when temperature dependencies (e.g. temperature- 
dependent viscosity) are included. 

Summing up, our numerical results are very promising. They raise hope that our proposed 
approach with the random aerodynamic drag force is capable of predicting the large elongations that 
are measured in industrial melt-blowing processes, presupposing an appropriate Cosserat model for 
the viscous, non-isothermal fiber jet. In addition, the computational effort seems to be manageable 
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Figure 4.9 

time scales tx — k/e and lx/v re / 



Adaptive time step choice At for (4.13) in comparison to turbulent 
k 3 ' 2 /(ev re j) with v re i = ||u — v|| experienced at 



the jet position r(t) of Fig. 4.6 



since the asymptotic globalization strategy for the sampling of the turbulent velocity fluctuations 
is linear in time and space discretization. 

Remark 11. Some concluding remarks on computational aspects: 



The adaptive time step control of the ODE-solver (cf. Remark 10) ensures the correct res 



olution of the turbulent scales since the chosen step size At is always clearly smaller than 
tx and lT/v re z, v re / = ||u — v||, cf. Figure 4-9 This implies the smooth numerical ap- 
proximation of the jet quantities, see e.g. the visualization of the velocity components in 
Figure \4-7\ 

The simulation of a fiber trajectory on [0, T), T = 10~ 3 [s] takes a CPU-time of approxi- 
mately 200 seconds on a 2.7 GHz Intel Core i5 processor. 

The computational effort of the sampling routine for the turbulent velocity fluctuations u' 
splits into initialization and continuous run. Whereas the costs for the initial generation 
of the set of random numbers are independent of the discretization and negligibly small 
(0.1 CPU-seconds for N — 50,), the costs for the continuous run are linear in the time 
discretization and add up to approximately 88% of the total costs for solving the random 
ODE system (4.13). On the first glance this seems to be incredibly much but the reason lies 



in the necessary processing of the underlying flow data (e.g. sorting, interpolation of flow 
data are required). So far, no further attention has been paid to the data processing that 
is done with standard MATLAB routines. But its performance will be optimized in future 
which promises a drastical speed-up. 



5. Conclusion and Outlook 

In a melt-blowing process liquid fiber jets arc spun due to turbulent air streams causing very high 
jet attenuation and final diameters of size smaller than a micrometer. So far, the understanding 
and design of the process suffered from a discrepancy between measurements and simulations; the 
computed final jet diameters were too thick by several orders of magnitude. This gap might be 
closed by considering the impact of the turbulent velocity fluctuations on the jets dynamic as 
we have demonstrated numerically in this paper. In correspondence to turbulence theory and on 
top of a k-e formulation we have modeled the turbulent velocity fluctuations as Gaussian random 
fields. Taking advantage of the special covariance/correlation structure we have proposed a fast and 
accurate sampling strategy whose effort is linear in the discretization and that makes the realization 
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possible. The numerical results are very convincing as they show already a qualitatively appropriate 
jet thinning in magnitude for a quite simple isothermal ODE-modcl for the jet dynamics. 

In future we intend to apply the developed random aerodynamic drag force to more sophis- 
ticated Cosserat models including also inner stresses and temperature dependencies in order to 
get quantitative estimates for melt-blowing. But this requires the robust numerical treatment of 
the PDE- models [2 El [6] . Especially, the handling of the expected huge elongations pose severe 
challenges on an efficient adaptive mesh refinement which is topic of recent research. 
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